#------------------------------------------------------------------------------
rm(list = ls())
library(LalRUtils)
libreq(
  data.table, zoo, tictoc, fixest, PanelMatch, patchwork,
  rio, magrittr, janitor, did, panelView, ggplot2, RPushbullet, ggiplot, 
  tidyverse, data.table, zoo, tictoc, fst, fixest, PanelMatch, patchwork,
  rio, magrittr, janitor, did, panelView, ggiplot, tictoc, binsreg, interflex
)

set.seed(42)
theme_set(lal_plot_theme())

notif = \(x) pbPost("note", x)

#------------------------------------------------------------------------------


#------------------------------------------------------------------------------
# Define Paths
#------------------------------------------------------------------------------

# R studio
setwd( dirname(rstudioapi::getActiveDocumentContext()$path) )
# R default : unccoment if you use default R
# setwd(getSrcDirectory(function(){})[1])
#------------------------------------------------------------------------------



#------------------------------------------------------------------------------
# Load Data
#------------------------------------------------------------------------------

vcf <- fread("vcf_data_complete.csv", sep = ",")
vcf[, "t" ] <- vcf$year - 1995
setnames( vcf, "d", "D")
vcf_data <- copy( vcf )
gfc <- fread("gfc_dta.csv" )


#------------------------------------------------------------------------------


#------------------------------------------------------------------------------
# Figure 4: Treatment Effects on Annual Deforestation as a 
#           Function of Ex-Ante Forest Cover Cutoffs
#------------------------------------------------------------------------------

library( purrr )
# Check the existance of vcf_data
# Filter data
if (exists("vcf_data")) {
  vcf = vcf_data[year >= 1995]
  rm(vcf_data)
}

# Generate subsample
# estimate with subsample using fixed effects and clusters
# Get results coeficients and standard errors
fitter = function(cutoff){
  dat = gfc[pref_bin >= cutoff]
  m = feols(def_ha ~ D | village[t] + styear, cluster = ~block, dat)
  tidy(m)[, 2:3]
}


fitter2 = function(cutoff){
  dat = gfc[pref_bin >= cutoff]
  m = feols(def_ha ~ D | village[t] + styear, cluster = ~block, dat)
  return(m)
}

models_gfc <- lapply( 1:10, fitter2 )

# Get results and add columns
tic()
cutoff_res = map_dfr(1:10, fitter) %>% setDT
toc()
cutoff_res$n = 1:10
colnames(cutoff_res)[1:2] = c('beta', 'se')

# GFC results
(
  rob_fit_gfc <- ggplot(cutoff_res, aes(n, beta)) +
    geom_point( size = 2 ) +
    geom_errorbar( aes(ymax = beta + 1.96 * se,
                       ymin = beta - 1.96 * se), 
                   alpha = 1, width = 0) +
    geom_hline( yintercept = 0 ) +
    scale_x_continuous( breaks = 1:10 ) +
    labs( title = 'GFC', 
          y = "Effect (Deforested Area)", 
          x = "Inclusion Threshold Decile of Forest Cover in 2000 ",
          caption = '' ) + 
    theme(axis.text = element_text( size = 22 ), 
          text = element_text( size = 23 ), 
          plot.title = element_text(size = 26), 
          legend.position = 'none')

)

# VCF Robustness
library(dplyr)
vcf[, pref_bin := ntile( cover_1990, 10 ) ]

# Function to subsample data
fitter = function(cut) {
  m = feols(forest_index ~ D | cellid[t] + styear, 
            data = vcf[pref_bin >= cut], 
            cluster = ~blk)
  tidy(m)[, 2:3]
}

fitter2 = function(cut) {
  m = feols(forest_index ~ D | cellid[t] + styear, 
            data = vcf[pref_bin >= cut], 
            cluster = ~blk)
  return(m)
}

# Get results and add columns
cutmods = map_dfr(1:10, fitter) %>% setDT
cutmods$n = 1:10
colnames(cutmods)[1:2] = c('beta', 'se')

models_vcf <- lapply( 1:10, fitter2 )

# VCF results
(
  rob_fit_vcf = ggplot( cutmods, aes( n, beta ) ) +
    geom_point( size = 2 ) +
    geom_errorbar(aes( ymax = beta + 1.96 * se,
                       ymin = beta - 1.96 * se), 
                  alpha = 1, width = 0) +
    geom_hline( yintercept = 0 ) +
    scale_x_continuous( breaks = 1:10 ) +
    scale_colour_brewer( palette = 'Set1' ) +
    labs( title = 'VCF',
          y = "Effect (Forest Index)", 
          x = "Inclusion Threshold Decile of Forest Cover in 1990",
          caption = '') + 
    theme(axis.text = element_text( size = 22 ), 
          text = element_text( size = 23 ), 
          plot.title = element_text(size = 26), 
          legend.position = 'none')
)

# Save figure
frob = ( rob_fit_vcf / rob_fit_gfc )
ggsave("main_figure4.pdf", frob, 
       device = cairo_pdf, 
       width = 8, height = 10 )

#------------------------------------------------------------------------------
